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ABSTRACT 

The intracluster medium (ICM) of galaxy clusters is a weakly collisional, high-beta plasma in which the 
transport of heat and momentum occurs primarily along magnetic-field lines. Anisotropic heat conduction 
allows convective instabilities to be driven by temperature gradients of either sign, the magnetothermal insta- 
bility (MTI) in the outskirts of non-isothermal clusters and the heat-flux buoyancy-driven instability (HBI) in 
their cooling cores. We employ the Athena magnetohydrodynamic code to investigate the nonlinear evolution 
of these instabilities, self-consistently including the effects of anisotropic viscosity (i.e. Braginskii pressure 
anisotropy), anisotropic conduction, and radiative cooling. We highlight the importance of the microscale in- 
stabilities (firehose, mirror) that inevitably accompany and regulate the pressure anisotropics generated by the 
HBI and MTI. We find that, in all but the innermost regions of cool-core clusters, anisotropic viscosity signif- 
icantly impairs the ability of the HBI to reorient magnetic -field lines orthogonal to the temperature gradient. 
Thus, while radio-mode feedback appears necessary in the central few tens of kpc, heat conduction may be 
capable of offsetting radiative losses throughout most of a cool core over a significant fraction of the Hubble 
time. Magnetically-aligned cold filaments are then able to form by local thermal instability. Viscous dissipation 
during the formation of a cold filament produces accompanying hot filaments, which can be searched for in 
deep Chandra observations of nearby cool-core clusters. In the case of the MTI, anisotropic viscosity main- 
tains the coherence of magnetic-field lines over larger distances than in the inviscid case, thereby providing a 
natural lower limit for the scale on which the field can fluctuate freely. In the nonlinear state, the magnetic 
field exhibits a folded structure in which the field-line curvature and field strength are anti-correlated. These 
results demonstrate that, if the HBI and MTI are relevant for shaping the properties of the ICM, one must 
self-consistently include anisotropic viscosity in order to obtain even qualitatively correct results. 
Subject headings: conduction — instabilities — magnetic fields — MHD — plasmas — galaxies: clusters: 
intracluster medium 

1. INTRODUCTION 

Clusters of galaxies are filled with hot and tenuous plasma, 
the intracluster medium (ICM), the detailed properties of 
which governs such important physics as heat and momen- 
tum transport, magnetogenesis, and thermodynamic stability. 
These properties are complicated by the fact that the ICM 
is only weakly collisional: while the particle mean free path 
Amfp is ^10-10'' times smaller than the thermal-pressure scale 
height H, it is neverthel ess ^lO'^-lO'^ times larger than the 
ion gyroradius r„\ (e.g. |Schekochihin & Cowley 2006 and 
references therein). As such, the material properties of the 
ICM are strongly anisotropic with respect to the magnetic- 
field direction, despite the fact that the strength of the intra- 
cluster magnetic field is relatively weak (^0.1-10 /iG, which 
constitutes only ^0.01-1% of the thermal energy; for a re- 
view, see Carilli & Taylor 2002j|. 

This anisotropy fundamentally changes the convective sta- 
bility properties of the ICM (Balbus 2000). Temperature 
gradients, rather than entropy gradients, become the dis- 
criminating quantities that determine stability, regardless of 
whether temperature increases ( Balbus, |2000| [2001 1 or de- 
creases (Quataert 2008) in the direction of gravity. As non- 
isothermal clusters generally exhibit both regions of positive 
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and negative temperature gradients, the entire ICM ought to 
be linearly unstable to convective motions. 

The temperature in the cores of non-isothermal clusters de- 
creases in the dkection of gra vity due to efficient radiative 
coohng at high densities (e .g. |Fabian|[T994[ |Piffaretti et al.| 
|2005[ jWdihnin et al.||2005"] l, and any alignment of conduct- 
ing magnetic-field lines and gravity th ere can lead to a heat- 
flux buoyancy-driven instability (HBI; Quataert'2008). Non- 
linear numerical simulations of the HBI have revealed that 
the instability acts in such a way as to quiescently shut itself 
off, gradually reorienting magnetic field to be perpendicular 
to the temperature gradient and thus stifling the heat flux that 



gave rise to the instability in the first place ( |Parrish & Quataert 
2008 , McCourt et al.|2 0lT). In the presence of radiative cool- 
ing, this field-line reorientation ultimatel y insulates the core, 
exacerbating the c ooling-flow problem (Parrish et al.||20()9{ 
[Bogdanovic et al.| |2009; Mikel lides et al.||2011| i u nlessleld 
lines are re-opened by sufficient turbulent stirring (Par rish e"t| 
al. 2010; Ruszkowski & Oh 2010; McCourt et al. 2011). 

Beyond the cooling radius, the temperature increases in the 
direction of gravity due to virialized gravitational infall. In 
this case, any m/ialignment of magnetic-field lines and grav- 
ity ca n lead to a magnetothermal instability (MTI; Balbus] 
|2000| |2001). In the presence of a sustained temperature gra- 
dient, the MTI leads to vigoro us subsonic turbulence and a 
radially biased magn etic field ( Parrish & Stone|[2005 2007 



McCourt et al.|20lT[ ). The former may provide up to 5-30% 
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of the pressure support beyond rsoo, reducing the observed 
Sunyaev-Z el'dovich signal and biasing X-ray mass estimates 
of clusters ( jParrish et al.|201 1} . The latter can lead to efficient 
radial heat transport, resulting in large-scale temperature pro- 
files flatter tha n those expecte d from structure formation cal- 
culations ( Panish et al.|2008^ . 

These studies of the HBI and MTI did not include an 
important feature of weakly collisional plasmas: changes 
in magnetic -field strength and/or density that occur on 
timescales much greater than the inverse of the cyclotron fre- 
quency result in pressure anisotropics (i.e. the gas pressure 
perpendicular and parallel to the local magnetic field become 
unequal). These pressure anisotropics manifest themselves as 
anisotropic viscous stresses, which target precisely those mo- 
tions originally responsible for the anisotropics themselves. 
By means of a linear stability analysis, which self-consistently 
accounted for the dynamical effect s of both anisotropic con- 
duction and viscosity, Kunz] ( 201 1 
I sub] I 



hereafter Kll) found that 



the HBI and MTI, when subject to the pressure anisotropics 
they induce, are qualitatively and quantitatively changed from 
what earlier studies had suggested. 

In brief, instabilities that depend upon the conver- 
gence/divergence of magnetic-field lines to generate unstable 
buoyant motions (the HBI) are suppressed over much of the 
wavenumber space, whereas those that are otherwise impeded 
by field-line convergence/divergence (the MTI) are strength- 
ened (Kl 1). This not only reduces HBI growth rates but also 
increases the wavelengths of the fastest-growing modes to 
^HBi/H ^ 0.2-1 (increasing outwards) for typical cool-core 
parameters. Taking into considerat ion the non-local nature of 
these modes. Latter & Kunz (20121 conjectured that the field- 
line insulation thought to be a nonlinear consequence of the 
HBI would be attenuated in all but the innermost ^20% of 
cluster cores. Perhaps not coincidentally, these regions tend 
to be dominated by strong radio-mode feedback from pow- 
erful central dominant galaxies (e.g. Burns 1990, Mittal et al. 



2009 Sun|2009[[Blanton et al |2010 ). The fastest-growing lin- 
ear MTI rnodes, on the other hand, escape the effects of pres- 
sure anisotropy by orientating their velocities perpendicular 
to the magnetic field. However, anisotropic viscosity couples 
Alfven and magnetosonic waves in such a way that damped 
slow-mode perturbations excite a buoyantly unstable Alfvenic 
response when the temperature increases in the direction of 
gravity. Consequently, many wavenumbers previously con- 
sidered MTI-stable or slow-growing are in fact maximally un- 
stable (see fig. 2 of Kll for an example). 

These changes raise a number of questions. Is the field-line 
insulation thought to be a nonlinear consequence of the HBI 
attenuated by anisotropic viscosity? If so, can anisotropic vis- 
cosity help conduction stave off a cooling catastrophe over as- 
trophysically relevant timescales? What role does anisotropic 
conduction and viscosity play in the generation of the cold fil- 
aments com monly observed in cluster cores (e.g. Lynds 1970; 
|Fabian eral.^20 08)? Does anisotropic viscosity significantly 
affect the ability of the HBI and MTI to amplify magnetic 
fields and drive turbulence? Are the resulting field strengths, 
magnetic-field topologies, and turbulent ve locities compati- 
ble with tho se observationally inferred (e.g. Schuecker et al.| 
[20041 [\^ & EnBHn 2005, Sanders et al.|2011| l? In this pa- 
per, we employ numerical simulations to understand these in- 
stabilities and the implications they have for the observable 
structure and evolution of a weakly collisional ICM. 

An outline of the paper is as follows. In Section |2] we 
present the basic equations describing a weakly collisional 



ICM, identifying the important dimensionless free parameters 
that govern the plasma physics. Section [3] describes our nu- 
merical approach and details our treatment of the microscale 
plasma instabilities that inevitably develop in our simulations. 
We then present our results concerning the non-radiative HBI 
(Section |4|, the radiative HBI (Section |5]l, and the MTI (Sec- 
tion |6]). Finally, in Section |7] we provide a brief summary of 
these results, a discussion ortheir astrophysical implications, 
a comparison with related work, and an outlook of what is 
required improve our understanding of the ICM. 

2. BASIC EQUATIONS 

The fundamental equations of motion may be written in 
conservative form as 



dp 
dt 



+ V-(pv) = 0, 



d(pv) 
dt 



+ V-(pvv+P*)=pg, 



dE 
Ik 



+ V ■ {Ev + P* -v + q) =pg-v-pC, 



dB 

Ik 



+ V-(vB-Bv) = 0, 



(1) 
(2) 
(3) 
(4) 



where p is the mass density, v is the velocity, g is the gravita- 
tional acceleration, B is the magnetic field. 



1 9 P B 
2^ 7-1 Stt 



(5) 



is the total (kinetic + internal + magnetic) energy density, and 
p is the gas pressure. The ratio of specific heats 7 = 5/3. 
We consider a hydrogenic plasma with equal ion and electron 
temperatures, 7] = Tg = T, and number densities, «i = n^, so 
that the mean mass per particle is nip/l. The thermal speed of 

the ions is then V(h = (pj pY^^ = {Ik^T /m^)^/^ . 

In Equations |2] and |3] we have introduced the total (gas + 
magnetic) pressure tensor 



g2 - 



(6) 



where p±^ (/^y) is the gas pressure perpendicular (parallel) to 

the magnetic field and b = B/B is the unit vector in the direc- 
tion of the magnetic field. The total gas pressure satisfies 



2 1 



(7) 



Differences between the perpendicular and parallel gas pres- 
sure arise from the conservation of the first and second adia- 
batic invariants for each particle on times cales much greate r 



than the inverse of the gyrofrequency, fig' (|Chew et al. 



19561. 



When the ion-ion collision frequency is much larger than 
the rates of change of all fields, an equation for the pressure 
anisotropy can be obtained by balancing its production by adi- 
abatic invariance with its relaxation via collisions: 



/7_L - P II =0.960 X '^-In^, 
" I'ii df 



(8) 



where pi is the ion gas pressure (e.g. |Catto & Simakov|2004[ ). 
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Defining the (ion) parallel viscous diffusivity, 



=0.960 X 



; 0.031 



ii 

2 Mi 



(9) 



Hi 



0.01 cm- 



and using Equations ([T]) and Q to replace the time derivatives 
of density and magnetic -fielcfstrength with velocity gradients, 
the pressure anisotropy (EquationlSll may be written 



P±-p\\ = 3pi^|| [bb--\ I :Vv 



(10) 



This pres sure an isotropy is the physical effect behind what is 
known as |Bragi nskii ( 1965) viscosity - the restriction of the 
viscous damping (to dominant order in the Larmor radius ex- 
pansion) to the motions and gradients parallel to the magnetic 
field. In an incompressible fluid, small-amplitude parallel- 
velocity fluctuations with parallel wavenumber are damped 
at a rate 



a;„isr = 3fcHi^ii 



(11) 



;0.61 



k\\H 

To" 



Mi 



0.01 cm- 



5/2 



2keV 



Motions that do not affect the magnetic-field strength to linear 
order (e.g. Alfven waves) are allowed at subviscous scales. In 
the weak-field regime, these motions take the form of plasma 
instabilities (see Section [J!2| ). 

The vast disparity between the gyro- and colUsion frequen- 
cies also implies that the heat flux q is anisotropic with respect 
to the magnetic field (Braginskii 1965) : 



9 = -pK||**- Vv^, 



(12) 



where 



Kii =1.581 X 



1.67 



1^ 

2 I^ee 
«i 

0.01 cm- 



(13) 



kBT 
2keV 



5/2 



kpc^ Myr ' 



is the (electron) thermal parallel diffusivity. 



: IkBT/lUe 



is the square of the electron thermal velo city, and i^ee is the 
electron-electron collision frequency (e.g. Catto & Simakov 
[2004) . Equation ( [T2| states that heat is transported along 
magnetic-field lines when there is a component of the temper- 
ature gradient aligned with the magnetic field. Field-aligned 
temperature fluctuations with parallel wavenumber fcy are dif- 
fused away at a rate 



t^cond = '^ll 



(14) 



i4.1 



k\\H 



0.01 cm- 



5/2 



Gyr-'. 
2 keV / ^ 



For future reference, we also define the (electron) parallel 
thermal conductivity x\\ = Pi^w/T T^^^. 

The ratio of the viscous and thermal diffusivities is known 
as the Prandtl number Pr, which is roughly constant: 



Pr; 



= 0.607 

K\\ Ai 



nil 



1/2 



: 0.02, 



(15) 



where Ag (Ai) is the electron (ion) Coulomb logarithm. This 
implies that viscous forces operate on a timescale that is a 
fixed number greater than the timescale on which conduction 
operates. In addition, there are two more important dimen- 
sionless parameters: the plasma beta. 



87r/7 



(16) 



and the inverse of the Knudsen number. 



far 

2keV 



Kn-' = 



H 



*mfp 



(17) 



1207 



10"^ cm s-2 



0.01 cm- 



2keV 



where H = v\/g is the thermal-pressure scale heightj^ The 
Knudsen number is a dimensionless measure of collisionality, 
and determines whether a fluid (rather than kinetic) descrip- 
tion may be used. Introducing the dynamical frequency 



Wdyn : 



;5.1 



1/2 



(18) 



10- 



2keV 



-1/2 



Gyr" 



the Knudsen number may equivalently be expressed as a ratio 
of frequencies: Kn = Wdyn/t'ii- In terms of Kn, the viscous and 
thermal parallel diffusivities are 



z/|l =0.48Knvth//, 



:24Knythi/, 



(19) 



(20) 



respectively. 

Finally, the last term in Equation [3] represents radiative 
losses. The cooling in the ICM is dominated by thermal 
Bremsstrahlung above temperatures ~'l keV, for which the 
radiative cooling rate (per unit volume) is 



pC ~ 10" 



»i 

0.01 cm- 



2keV 



1/2 



ergs cm ^ s ' 



(21) 



(j Rybicki & Lightman 19791. For an isobaric perturbation, 
Equations and (|21[) imply a cooUng frequency 



Wcool : 



2 d{pC/p) 



5 d\nT 



(22) 



-0.3 



0.01 cm- 



kuT 
2keV 



-1/2 



Gyr"'. 



The fact that the cooling frequency is negative indicates that 
isobaric thermal instability is possible (Field 1965). Conduc- 
tion suppresses this instability (at least along field lines) for 

^ In determining tlie numerical value of Kn, we have assumed force bal- 
ance between thermal pressure and gravity - an assumption that will hold as 
an initial condition in all of our simulations. 
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parallel wavelengths smaller than the Field length. 



AF = 27r - 



2 Kl] 



:300 



5 Wcool 



1/2 



(23) 



1 / k^T 



3/2 



kpc. 



.0.01cm-3y V2keVy 
We employ radiative cooling in two of our HBI simulations. 

3. NUMERICAL APPROACH 

3.1. Integration Scheme 

We integra te equati ons (fTll-(j4|i using the conservative MHD 
code Athena ([Stone et al. 2D08 ). Details concerning the MHD 
algorithms may be found in Gardiner & Stone (2005 2008). 
The directionally unsplit corner transport upwind (CTU) inte- 
gration method and the Roe Ri emann solver are used in all 
of our simulations. Following Sharma & Hammett (2007i 



and Dong & Stone (2009 ), respectively, anisotropic conduc- 
tion and Braginskii viscosity are implemented via operator 
splitting using slope limiters on the transverse heat and vis- 
cous fluxes to ensure stability. The conduction algorithm is 
sub-cycled with respect to the main integrator with a time 
step Afcond = C{Axf/[2d('y- 1)k||], where C w 0.5 is the 
Courant number and d is the number of spatial dimensions be- 
ing solved. In order to prevent impulsive driving due to abrupt 
changes in pressure, we restrict the sub-cycling routine to take 
no more than 10 sub-cycles per global timestep. 

When radiative cooling is included in our study of the HBI 
(see Section |5|, we employ the exact integration scheme de- 
tailed in Townsend (2009 ). The cooling source term is added 
to the reconstruction and interface-state correction steps in the 
CTU integrator, so that the cooling is fully second-order ac- 
curate. In order to prevent the for mation of an unresolved 
cold phase, we follow jSharma et ar| ( j2010^ and McCourt et al.| 
(^012) in adopting a temperature floor Tfloor = 7o/20, where Tq 
is the initial minimum temperature in our model atmosphere. 
This temperature floor is based on the reasonable assump- 
tion that, once a thermally unstable fluid element cools below 
Tfloor, it is unlikely to enter back into the hot phase. Because 
our focus is on the evolution of the HBI, subject to radiative 
cooling and Braginskii viscosity, and not on the detailed na- 
ture of multiphase gas in a thermally unstable ICM, this sim- 
pUfication should not significantly affect our results. 

3.2. Microscale Instabilities 
When the pressure anisotropy violates the inequalities 



B2 

47r 



If_ 
8^' 



(24) 



rapidly growing microscale instabilities (firehose and mirror, 
respectively) are trigge red and the Braginski i-MHD equations 
become ill-posed (see |Schekochihin et ar][2005 and refer- 
ences therein). Without finite Larmor radius effects taken into 
account, the fastest growing microscale modes formally oc- 
cur at infinitely small scales, which in practice translates to 
scales near the grid where the microscale instabilities may 
be unresolved. Exactly what to do in this situation is not 
obvious and is currently under investigation (A. Schekochi- 
hin, private communication). In the mean time, because the 
pressure anisotropy controls the rate of viscous dissipation - 
which in turn affects the large-scale dynamics - some mea- 
sures must be taken in order to capture the microscale in- 



fluence on the pressure anisotropy, particularly in the weak- 
field regime. In this paper we choose two approaches, both of 
which are supported by strong evidence in the solar wind and 
magnetosheath (plasmas in many ways similar to the ICM) 
that microinstabilities isotropize the plasma to marginally- 
stable levels (e.g. Stevens & Kasper_2007 , Bale et al. 2009) . 

Our first approach is based upon on the theory that, once 
triggered, microscale fluctuations grow in such a way as to 
compensate on average the "excess" pressure anisotropy gen- 
erated by the large-s cale motions, thus maintaining ma rginal 
stability ( Schekochi hinet al.| j2008UR osin et al.| [20TT). We 
simply allow the microscale mstabilities to self-consistently 
develop over the course of our simulations and to naturally 
regulate the pressure anisotropy, with the expectation that nu- 
merical viscosity will prevent the relatively small structures 
from getting out of hand. Our simulations are deliberately 
chosen with high enough resolution to not only resolve the 
HBI and MTI, but also to ensure healthy time- and lengthscale 
separations between the HBI/MTI and the firehose/mirror in- 
stabilities. As a result, the firehose fluctuations we resolve in 
our simulations grow fast enough to rapidly enforce marginal 
stability and self-consistently provide a hard-wall limiter on 
negative pressure anisotropics. Unfortunately, the same can- 
not be said for the mirror instability. The Braginskii version of 
the mirror instability grows substantially slower than the ki- 
netic mirror instability, having a growth rate smaller than the 
parallel rate-of-strain of the viscous scale eddies. While neg- 
ative pressure anisotropics are efficiently regulated, positive 
pressure anisotropics may not be. 

It is important to note that, because we are not able to 
simultaneously resolve the ion Larmor radius and thermal- 
pressure scale height, which requires >13 orders of magni- 
tude in scale separation, the microscale instabilities triggered 
throughout the course of our simulations do not grow as fast 
as they would otherwise grow in nature. In our Braginskii- 
MHD simulations the maximum growth rate of the firehose 
instability, fe||.n,axVth|A + 2//3p/2 A^Wdyn|A + 2//3|'/2, occurs 
at k\\H ^ N, where A = (p±-p\\)/p is the fractional pressure 
anisotropy and is the number of grid zones per thermal- 
pressure scale height. By contrast, a kinetic calculation in- 
cluding FLR effects reveals that the parallel firehose actu- 
ally has a maximum growth rate ^rigi|A + 2//3| occurring at 

A:||rg_i ^ |A + 2//3|'/^, spreading to larger scales as the pres- 
sure anisotropy approaches marginality. Even for our highest- 
resolution simulation (A^ = 512), we are underestimating the 
maximum growth rate of the firehose instability by a factor 
--10"|A + 2/,3|'/^. This is one example of the fact that the 
nonlinear saturation of microscale fluctuations and the con- 
sequent regulation of pressure anisotropy occurs in our sim- 
ulations on a timescale much longer than it would in na- 
ture, where microscale fluctuations grow to SB /B ^ I on a 
timescale comparable to the turnover time of the turbulent 
motions. A potentially serious consequence of not resolv- 
ing the microscale instabilities at their natural scales is that 
we may be overestimating the conductivity and viscosity of 
the plas ma by a factor ^Amfp/ yg.i ^ 10'°-10" (see final para- 
graph of |Schekochihin et al.|2008| ). 

Our secon d approach is motivated by the work of [Sharma] 
et al. ( 2006"), who numerically investigated the nonlinear evo- 



lution of the collisionless magnetorotational instability and 
accounted for the effects of microscale instabilities by artifi- 
cially limiting the pressure anisotropy to lie within the bounds 
given by Equation ([24|. The computational advantage of 
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this approach is that microscale fluctuations are never trig- 
gered during the simulation. This closure rests on the follow- 



instabilities are triggered: 



ing plasma-physical rationale (e.g. Schekochihin & Cowley 
|2006 ). Once these thresholds are crossed, microscale insta- 
bilities will produce a fluctuation "foam" off of which parti- 
cles may pitch-angle scatter, break adiabatic invariance on the 
extremely short cyclotron timescale, and thereby isotropize 
the pressure (provide d such fluctua t ions c an penetrate down 
to the ion gyroscale). Sharma et al. ( 2006| l modeled this pro- 
cess by a large effective collisionality, which was activated in 
regions where and when the microscale stability boundaries 
were sufficiently exceeded, its magnitude being proportional 
to the product of a large frequency f p ^ 1 /A? and the pres- 
sure anisotropy excess. This effectively raises the Reynolds 
number of the plasma and makes it more coUisional. 

Note that when and where this occurs the pressure 
anisotropy is no longer connected to the large-scale turbulent 
stretching of the magnetic field that gave rise to the pressure 
anisotropy in the first place, nor does this approach capture 
the effects of the microscale contribution to the total rate-of- 
strain of the plasma. Moreover, the heating associated with re- 
laxation of the pressure anisotropy is not correctly captured, 
as the assumed rapid pitch-angle scattering and consequent 
pressure isotropization has no associated heating term in our 
energy equation. Indeed, d eciding exactly wha t to do with 
this energy is not trivial (see Sharma et al.|2007 1 



In summary, our two approaches amount to two different 
interpretations of Equation ^ in the presence of microscale 
instabilities, with antithetical implications for the viscous dis- 
sipation of macroscale motions. In the first approach, the 
collision frequency of the plasma remains constant while the 
microscale instabilities modify on the average the (parallel) 
rate-of-strain so as to offset the pressure anisotropy caused 
by the changing macroscale fields (i.e. dlnZ?/df oc va/ P). At 
the macroscales, the plasma behaves as though it were more 
viscous. In the second approach, the microscale instabilities 
break adiabatic invariance, effectively increasing the colli- 
sion frequency { va — > i^p ^ Wdyn/?) and returning the pressure 
anisotropy to marginally stable values. At the macroscales, 
the plasma behaves as though it were less viscous. The im- 
portant question of which of these interpretations is correct 
boils down to the (unanswered) question of whether or not 
such microscale fluctuations can reach the ion Larmor radius 
in a driven, initially Maxwellian system. 



3.3. Choice of Dimensionless Plasma Parameters 

Our choice of dimensionless parameters is motivated by 
considerations of both the physical conditions in actual galaxy 
clusters and the numerical constraints related to the above 
pressure-anisotropy concerns. Since short-wavelength modes 
with k\\H > Z?'/^ are stabilized by magnetic tension (unless 
A + 2//3 < 0), one would ideally like to construct simulations 
with relatively large (3 so that a healthy spectrum of HBI and 
MTI modes may grow unabated (at least in their linear phase). 
However, such /? not only are much larger than the observa- 
tionally estimated ICM /3 ^ 10^-10^ (for a review, see Cariliil 
|& Tay lor 2002 ), but also place steep constraints on how long 
the HBI and MTI can be evolved without our simulation re- 
sults being plagued by the aforementioned microphysical un- 
certainties. One can estimate from Equations (jSll and ( 24 1 how 
large SBiJB can grow in the linear phase before microscale 



^ l_ 

B /3f7HBi ^ /3Kn' 



(25) 



In cluster cores Kn ^ 10~^-10~^, increasing outwards, and a 
choice of very large /? implies that one cannot go far beyond 
the linear regime without running the risk that the microphysi- 
cal closures we have employed significantly influence the sub- 
sequent large-scale dynamics. We therefore choose (3 ^ IQ*- 
10^ in our simulations, which puts them in a regime in which 
Braginskii viscosity is more important than magnetic tension 
and in which the instability can develop SB^^/B ^ a few per- 
cent before firehose and mirror instabilities are triggered. In 
the outer regions of galaxy clusters, Kn ^ 10~^-10~' and it is 
almost trivial to violate Equation ( [25] l. 

All this being said, we believe that our simulation results 
represent a step forward in our understanding of convective 
instability and thermal conduction in the ICM. In lieu of a full 
kinetic simulation that can resolve >13 orders of magnitude 
in spatial and temporal scale or a sub-grid model that can cor- 
rectly capture the complex interplay between the micro- and 
macroscales, this seems to be the best we can do at this stage. 

4. NON-RADIATIVE HBI 

4. 1 . Background Equilibrium and Initial Perturbations 

We consider a non-radiative, plane-parallel plasma strat- 
ified in both density and temperature in the presence of 
a uniform gravitational acceleration in the vertical direc- 
tion, g = -gz. The plasma is threaded by a uniform back- 
ground magnetic field oriented along z and is assumed ini- 
tially Maxwellian so that p± = = p in the background state. 
Force balance then implies 



dp 
dz 



= -gP- 



There is a heat flux in the background state given by 

dr. 

9 = -X||^^- 



(26) 



(27) 



In order to preserve thermal equilibrium, the background heat 
flux must be divergence-free: 



dzV^"S)='- 



(28) 



Enforcing T = To at z = Q and T = Tz at z = Z, Equation ( 28 1 
may be integrated to yield the temperature profile 



T(z) = Tc 



2/7 



(29) 



where ( = {Tz/Tq)' \ measures the magnitude of the steady 
heat flux through the atmosphere {q oc Q. Combining this 
result with Equation ([26]l determines the pressure profile 



/?(z) = /?oexp 



7 Z 
"5C^ 



1 + 



5/7 



1 



(30) 



where Hi 







''th.O/ 



is the thermal-pressure scale height and 
P() is the thermal pressure, both evaluated at z = 0. Note that 
dlnT/dz is largest at z = 0, and so HBI modes will naturally 
grow fastest at small z. 
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Table 1 

Parameters Used in HBI and MTI Simulations 



Run 


d 


Box Size 


Resolution 


Field Configuration 


Kn;:' 




Notes Section 


H2dBrag 


2 


H(\ X 2H{\ 


512x 1024 


vertical; /3o = 


105 


1500 


Braginskii viscosity 


4.3 




H2dlsoP 


2 


HqxIHo 


512x1024 


vertical; /3o = 


105 


1500 


isotropic pressure 


4.3 




H2dBLim 


2 


HoXlHo 


512x1024 


vertical; /3o = 


105 


1500 


artificially limited Braginskii viscosity 


4.3 




H3dBrag 


3 


HoxHoXlHo 


128x128x256 


vertical; /3o = 


10= 


1500 


Braginskii viscosity 


4.4 




H2dBRad 


2 


HoXlHo 


512x1024 


vertical; /3o = 


105 


1500 


Braginskii viscosity and thermal Bremsstrahlung 


5 




H2dIRad 


2 


UqxIHo 


512x1024 


vertical; = 


10= 


1500 


isotropic pressure and thermal Bremsstrahlung 


5 




M2dBrag 


2 


HqxIHo 


512x1024 


horizontal; /3o 


= 105 


200 


Braginskii viscosity 


6.3 




M2dlsoP 


2 


HoXlHo 


512x1024 


horizontal; /3o 


= 10= 


200 


isotropic pressure 


6.3 




M2dBLim 


2 


HoXlHo 


512x1024 


horizontal; /3o 


= 105 


200 


artificially limited Braginskii viscosity 


6.3 




M3dBrag 


3 


HqxHqx2Ho 


128x128x256 


horizontal; /3o 


= 10^ 


200 


Braginskii viscosity 


6.4 




M3dIsoP 


3 


HqxHqxIHo 


128x128x256 


horizontal; /3o 


= 10^ 


200 


isotropic pressure 


6.4 





This equilibrium is characterized by two dimensionless free 
parameters: C,, which is ~' 10-50 in cool-core clusters, and 



Ho 



(31) 



V250kpc 



2keV / 



which is a measure of the height of the atmosphere. We 
choose Q = 2.0 and Tz/To = 2.5 (i.e. C — 23.7), which implies 
Pzl Po — 0.14 by Equation (30 1. These numbers are charac- 
teristic of the cool-core cluster A 1795 with k^To, k, 2.5 keV, 



k^iTz ~ 6.3 keV, and Z w 250 kpc ( |Ettori et al.|2002 1. 

We apply Gaussian-random velocity perturbations to our 
background equilibrium, having a flat spatial power spectrum 
and a standard deviation of 10~^Vth.o- Such perturbations are 
sufficiently subsonic to ensure linear evolution from the out- 
set. These initial conditions are not representative of those 
conditions found in actual clusters, in which galaxy motions, 
major and/or minor mergers, and feedback from active galac- 
tic nuclei (AGN) stir the plasma. 

4.2. Numerical Setup and Boundary Conditions 

The equations are put in dimensionless form by choosing 
units natural to the problem. The units of velocity [v], density 
[p], and magnetic -field strength [B] are, respectively, the ini- 
tial values of the thermal speed v'th.o, density pQ, and magnetic 
field Bo at the bottom of the box (z = 0). The unit of length [I] 
is //(), so that the implied unit of time [t] is //o/vth,o, the initial 
sound-crossing time across a thermal-pressure scale height. 
Since pressure balance applies in the equilibrium state, g=\ 
in these units. The units of length and time have the scalings 



{£]: 



and 



[f] ~ 196 



124 f 

\\Q-^ cm s-2y 

(lO-8 cm s-z) ( 



kuTj) 
2keV 



1/2 



kpc (32) 



2kev) ^y-"' 



(33) 



respectively. Aside from the free parameters C and Q associ- 
ated with the background equilibrium, the linear evolution of 
the HBI depends only upon the plasma beta and the Knudsen 
number We choose /3o = 10^ and Kjiy' = 1500, so that Bq ~ 
0.016, v\\ ~ 3.2 X 10-"* T^l^p-\ and K|| ~ 1.6 x lO'^ T^/^p-i 
in dimensionless units. These imply initial values of v\\ ~ 



0.023 and K|| ~ 1.1 at z = Z; the latter causes our simulations 
to be rather expensive. 

Linear analysis of the HBI with Braginskii viscosity has 
shown that the only modes to evade strong suppression are 
confined to a thin band in wavenumber space in which con- 
duction is fast but viscous damping is small: Wcond ^ Wdyn ^ 

i^visc, or, using the definitions |TT| and |l4| l, k\\H < a/I^ < 
3A;||// (Kl 1). Taking the initial temperature and pressure pro- 
files. Equations (29 1 and (|30]) respectively, and using equa- 
tions (49) and (50) from kTI, we expect the fastest growth 
at small z on parallel wavelengths A|| w O.IHq. However, 
the lower coUisionaUty at larger z shifts the fastest-growing 
modes to w avelengths comparable to the thermal-pressure 
scale height ( Latter & Kunz|[2"012| l. In order to capture this 
global behavior, we choose box sizes HqxIHo (in 2d) and 
HoxHoxlHo (in 3d), so that L, = Z = 2Hq. We have also 
run local simulations with the size of the simulation domain 
much smaller than the ther mal-pressure scale hei ght similar 
to those presented in §4.1 of'McCourt et al.| ( |201 1| and found 
that only extremely slow-growing HBI niodes fit into the box, 
in agreement with linear theory. Our 2d runs have a resolu- 
tion of 512x1024. Our 3d run necessarily has lower reso- 
lution (128x 128x256) due to the stiff numerical constraints 
imposed by heat and momentum diffusion in three dimensions 
(a single 3d run at this resolution requires ^40,0 00 CPU-hrs). 
T he b oundary conditions are the same as in [McCourt et| 



al. (2011). The temperatures at the upper and lower bound- 



aries of our computational domain are fixed for all times, 
T{z = 0,f) = 7b and Tiz = Z,t) = Tz, while the pressure is 
extrapolated into the upper and lower ghost zones in such 
a way as to ensure hydrostatic equilibrium at those bound- 
aries. These choices are motivated by the observation that 
many galaxy clusters in the local universe are observed to 
have non-negli gible temperatur e gradients (e.g. Piffaretti et 
|al. |j j2005] iVikhlinin et al.|[20a5] ). The magnetic field is con- 
strained to cross the upper and lower boundaries normally, al- 
though its strength there is allowed to adjust according to the 
local dynamics. Periodic boundary conditions are imposed in 
the horizontal direction(s). 

Here we present results from four non-radiative HBI sim- 
ulations. H2dIsoP is a 2d simulation with isotropic pressure 
and serves as a reference run, enabling us to draw conclu- 
sions about the effects of Braginskii viscosity. H2dBrag is a 
2d simulation with Braginskii viscosity and H2dBLim is a 2d 
simulation in which the pressure anisotropy is artificially lim 



ited using the Sharma et al. closure described in Section 3.2 
H3dBrag is a 3d simulation with Braginskii viscosity. The 
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Figure 1. Temporal evolution of the box-averaged horizontal (x) and ver- 
tical (z) kinetic and inagnetic energy densities in runs H2dBrag (red lines), 
H2dIsoP (blue lines), and H2dBLim (purple lines). The units of energy den- 
sity and time are, respectively, pov^j, q and //o/i'ih,o (see Equationjssj. 

parameters in these simulations are summarized in Table [T] 

4.3. 2d Simulations 

In Figure[T]we present the evolution of the box-averaged ki- 
netic and magnetic energy densities. Runs in which pressure 
anisotropics are allowed to develop (red and purple lines) ini- 
tially exhibit a growth rate «0.6 times smaller than that of the 
run without Braginskii viscosity (blue line), in agreement with 
predictions from linear theory (Kll). A comparison between 
the growth rates of individual Fourier modes in run H2dBrag 
and those predicted by a quasi-global linear stability analy- 
sis of our model atmosphere shows e xcellent agreement al l 
the way to k,H - 3000 (see fig. 5 of [Latter & Kunzl[20l2l l. 
The growth rate in run H2dBLim (purple line) departs from 
that of runH2dBrag (red line) once the hard-wall pressure- 
anisotropy limiters become active and regulate the pressure 
anisotropy. The kinetic energy thereafter grows similarly to 
the run without Braginskii viscosity. All three runs reach an 
approximately saturated kinetic energy density corresponding 
to a box-averaged Mach number of a few percent (although 
local velocities can range up to «50% of the local thermal 
speed). The total magnetic energy increases by a factor of 
^10 over the course of the runs, similar to the increase in the 
total kinetic energy during the non-exponential phase of evo- 
lution. Horizontal and vertical energies reach approximate 
equipartition by the end of the simulation. 

While all three runs appear similar, box-averaged quantities 
can be deceiving. The overall spatial and temporal evolution 
of the atmosphere during each of our 2d runs is shown in Fig- 



urel2] The temperature (color) and magnetic-field lines (black 
lines) are displayed in each of the six frames, which show the 
atmosphere at the different times t = 10, 15, 20, 30, 40, and 
50 (in units of Ho/vth.o; see Equation 33 1. 

In each of the runs, the HBI develops first at low altitude 
where the temperature gradient is largest, eventually progress- 
ing to higher altitudes where the temperature gradient is shal- 
lower. Runs with Braginskii viscosity (top and bottom rows) 
show a significant delay in the development of the HBI, par- 
ticularly at higher altitudes. There are three reasons for this 
difference. First, early viscous damping of the seed velocity 
perturbations causes the instability to grow from smaller am- 
plitudes than in the run with isotropic pressure. Second, vis- 
cous damping of motions that compress and rarefy the mag- 
netic field lines results in a reduced HBI growth rate. Third, 
the increasing importance of viscous damping with height (re- 
call i^ii cx r^/^p"') shifts the HBI to successfully larger paral- 
lel wavelengths, which eventually become comparable to the 
local scale height and lead to secular (rather than purely expo- 
nential) growth due to non-local effects. As a result, there is 
substantial difference in the saturated-state temperature pro- 
file and the topology of the magnetic-field lines for z > O.IL . 

In run H2dBrag, regions of negative pressure anisotropy 
(i.e. decreasing magnetic -field strength) produce firehose in- 
stabilities near the grid scale when the local value of {p±^ - 
Pj[ /An) < 0. The firehose fluctuations grow exponentially 
until they compensate for the excess pressure anisotropy, after 
which they grow secularly. Once the local pressure anisotropy 
is regulated to be ~-B-/47r, the firehose turbulence moves to 
longer wavelengths. This can be seen clearly by comparing 
the structure of the firehose modes near the top of the f = 50 
panel with those of the f = 20 panel. In run H2dBLim, such 
firehose fluctuations do not exist by construction: the Bragin- 
skii pressure anisotropy is limited by hand to lie within the 
microscale stability boundaries. Rather than produce firehose 
instabilities, regions of large negative pressure anisotropy ef- 
fectively eliminate the magnetic tension. As a result, sharp 
folds in magnetic fields are allowed to develop in regions of 
decreasing magnetic-field strength. Note that there is much 
less numerical reconnection in runs H2dBrag and H2dBLim 
than in run H2dIsoP, since small-scale motions that change 
the magnetic-field strength are viscously damped. 

The horizontally averaged magnetic-field angle as a func- 
tion of height, calculated as (/7?)j^, is shown in Figure b\ at 
the same times as in Figure [2] In all three cases (HldSrag: 
red line; H2dIsoP: blue line; H2dBLim: purple line), the HBI 
grows by reorienting the magnetic field to be more and more 
horizontal and {b^)x generally decreases in time at all heights. 
However, there are significant differences between each of the 
runs. As was evidenced in Figure [2] Braginskii viscosity im- 
pedes appreciable field-line reorientation for z > 0. ILj, where 
the viscous and dynamical frequencies become comparable. 
At these heights, the large (parallel) wavenumbers required to 
keep the HBI in action as the magnetic field becomes more 
and more horizontal are strongly suppressed by the pressure 
anisotropy they generate. For example, at f = 30 (~6.6 Gyr 
for g = 10"^ cm s"^ and k^To = 2.5 keV) the magnetic-field an- 
gle in the Braginskii and non-Braginskii runs are similar only 
in the innermost w20% of the core. Beyond this height, the 
relatively straight field lines in run H2dBrag allow heat con- 
duction to remain active at a rate comparable to the field-free 
Spitzer value. The magnetic-field angle in run H2dBLim, in 
which the pressure anisotropy was artificially limited so as to 
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Figure 2. Spatial and temporal evolution of the HBI with Braginskii viscosity (top row), without Braginskii viscosity (middle row), and with limited Braginskii 
viscosity (bottom row) in two dim ensions. The temperature (color) and magnetic-field lines (black lines) are shown at times ? = 10, 15, 20, 30, 40, and 50 (in 
units of //o/i'ih.o; see Equation |33| . The computational domain has dimensions L^xL- = 1 x2 (in units of Hq; see Equation |32| . We have suppressed the imaging 
of temperatures beyond the fixea color-bar limits. 



prevent microscale instabilities from growing, is intermediate 
between that of run H2dBrag and H2dIsoP. This is because 
the limiters restrict how effective Braginskii viscosity can be 
at suppressing the HBI. 

In Figure |4] we plot the box-averaged pressure anisotropy 
as a function of time in run H2dBrag. The solid (dashed) 
black line denotes an average positive (negative) pressure 
anisotropy, while the dotted line traces the box-averaged value 
of B^/A tt (a quantitative measure of microscale stability; see 
Section [372] ). Initially, the pressure anisotropy grows expo- 
nentially because pj_-p\\ oc SB^^/Bq. During this phase, there 
are more regions of decreasing field strength (^B|| < 0) than 
increasing field strength {SB\\ > 0) and so the box-averaged 
pressure anisotropy is negative. This is because regions with 
(5B|| < correspond to downward displacements, which pre- 
dominate by taking advantage of the steeper temperature pro- 



file at smaller z (note that d^T/dz^ < 0). Once these regions 
of negative pressure anisotropy satisfy p±-p\\ < -B\/Ati, 
rapidly growing firehose fluctuations efficiently reduce the 
pressure anisotropy to marginal stability. This is why the 
dashed line in Figure [4] never appreciably crosses the dot- 
ted line. Once the HBI settles into its nonlinear phase, the 
box-averaged pressure anisotropy is positive since, in gen- 
eral, p±-p\\ oc dlnB/df > 0. Because the mirror instability is 
not accurately captured in our Braginskii-MHD simulations, 
it is not as efficient at regulating the pressure anisotropy as it 
would be in nature. Nevertheless, the box-averaged pressure 
anisotropy always stays within a factor of a few of B^/4-tt, due 
to efficient firehose and not-so-efficient mirror regularization. 

We note here that the evolution of the kinetic energy shown 
in our Figure [T| is qualitatively differ ent than that presented 
in figure 3 of McCourt et al. ( 2011| l. Those authors found 
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Figure 3. Spatial and temporal evolution of the horizontally averaged 
magnetic-field angle at each height in runs H2dBrag (red lines), H2dIsoP 
(blue lines), and H2dBLim (p urpl e lines). The units of lengt h an d time are, 
respectively, Hq (see Equation|32l and //o/^'ih.o (see Equation|33l. 
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Figure 4. Temporal evolution of the box-averaged pressure anisotropy in run 
H2dBrag. Positive (negative) pressure anisotropics are denoted by a black 
solid (dashed) line. The dotted line traces the box-average d va lue of B~ /Att, 
a quantitative measure of microscale stability (see Section [X2) . The units of 
pressure and time are, respectively, poif), o ^'^^ / I'th.o (see Equation |33| . 

that the energy in the vertical motion was in the form of sta- 
ble oscillations that decay nonlinearly, whereas the horizon- 
tal kinetic energy persisted once the magnetic field became 
predominantly horizontal. While it is not surprising that our 
Braginskii-HBI simulations do not show this tendency, it is 
rather intriguing that our non-Braginskii-HBI simulations do 
not either We attribute this to their choice of = 10'^. At 
such small magnetic -field strengths, the magnetic field ex- 
erts essentially no dynamical feedback upon the gas motions, 
even as the magnetic field acquires a sharp folded structure. 
A more careful and dedicated study of this difference will be 
presented in Avara et al. (2012, in preparation). 

4.4. 3d Simulation 

Figure |3]presents the evolution of the box-averaged kinetic 
and magnetic energy densities from run H3dBrag. There are 
a few differences between the behavior shown in this figure 
and that shown in its 2d analog (Figure [TJ- Some of these 
are genuine differences due to the increased dimensionality 
of the simulation, while others are attributable to the factor of 
4 difference in resolution. To determine which differences are 
genuine and which are not, we have run a 128x256 version 
of run H2dBrag that is similar to run H3dBrag in every way 
except dimensionality. Comparing these two runs with the 
original 512x 1024 run H2dBrag, we are able to conclude that 
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Figure 5. Temporal evolution of the box-averaged kinetic and magnetic en- 
ergies in run H3dBrag. The units of en ergy density and time are, respectively, 
Po^'ii.o Ho/v(h.o (see EquationpsJ. 
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Figure 6. Spatial and temporal evolution of the horizontally averaged 
magnetic-field angle at each height in r un H 3dBrag. The units of lengt h an d 
time are, respectively, //q (see Equation|32| and //o/i'th.o (see Equation|33|. 



run H3dBrag is not fully converged (the energies in both the 
2d and 3d runs are generally too small by a factor of a few). 
While this is a lesson in itself - that properly simulating the 
HBI under actual cluster conditions requires more than 128 
zones per thermal-pressure scale height - we believe there is 
much to be gained from nevertheless presenting our 3d results. 

Of those differences that are clearly attributable to the in- 
creased dimensionality, one is a marked deficit of energy 
in the horizontal components of the velocity and magnetic 
field. This is because small-wavelength perturbations whose 
wavevectors have a component perpendicular to both grav- 
ity and the magnetic field behave like modified Alfven waves 
that are only slowly growing or decaying (depending on their 
exact wavevector orientation; see §4.1.2 of Kll). In other 
words, as the magnetic field becomes on the average more 
horizontal, the extra degree of freedom allows Braginskii vis- 
cosity to reorient these perturbations so as to minimize field- 
line compressions and rarefactions. 

Another difference is in the horizontally averaged 
magnetic-field angle as a function of height and time, shown 
in Figure |6] While the curves for f < 40 are very similar to 
those shown as red lines in Figure [3] (run H2dBrag), there is 
one noticeable difference: the magnetic-field angle, {b^)x.y, is 
smaller for z < 0. 1 in the 3d run. We attribute this to inter- 
change motions, which allow horizontally inclined magnetic 
field lines to slip past one another (a similar effe ct occurs in 3d 
simulations of the Rayleigh-Taylor instability; |Stone & Gar-| 
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diner 2007 1. It is in this region where the field lines have been 
considerably reoriented, since the effect of Braginskii viscos- 
ity is greatly reduced there (recall cx r^/^p"'). Elsewhere 
in the atmosphere, these interchange motions do not occur 
as readily, since the magnetic field does not become predom- 
inantly horizontal due to effective parallel-viscous suppres- 
sion of short-wavelength HBl modes. While the f = 50 panel 
in Figure |6] appears different than the red line in the t = 50 
panel of Figure with a relative decrease in {b^)x,y at all 
heights, we can safely attribute this to insufficient resolution: 
our 128x256 2d test simulation shows a similar decrease. 

These properties can be also be seen in Figure |7] which 
shows the magnetic-field lines (color-coded accordmg to the 
local field strength) at times t = 30 and 50. At f = 30 the HBI 
has yet to significantly affect the magnetic field for z^Hq. 
For O.lZ/o ^ z ^ O.6//0, the field lines have been strongly 
compressed in some regions and strongly rarefied in others. 
The flow of heat across these heights occurs along magnetic 
sheaths (or filaments), inside of which j3 ^ 10^. The resulting 
tension in these strong-field regions is responsible for straight- 
ening these field lines out. For z < O.I//0, the HBl has suc- 
cessfully reoriented the field lines to be predominantly hori- 
zontal. By f = 50, the HBl is active to various extents through- 
out the entire atmosphere. The magnetic filaments have be- 
come longer and more prominent, while the field lines below 
z ~ O.I//0 remain horizontally inclined. These magnetic bun- 
dles were also seen in our 2d simulations, but are more pro- 
nounced here in 3d. The reason is that it is more efficient to 
gather field lines together in 3d than in 2d. 

5. RADIATIVE HBI 

We have also run two 2d HBI simulations including radia- 
tive cooling, one with Braginskii viscosity (H2dBRad) and 
one with isotropic pressure (H2dIRad). The parameters used 
are summarized in Table [T] 

Our numerical setup is the same as that detailed in Section 
[4j aside from two important differences. First, when radiative 
cooling is considered, thermodynamic equilibrium requires 



t = 30 



d 

d~z 



XII 



dT 
dz 



= pCoc p T 



27,1/2 



(34) 



In contrast with Equation ( |28) , this equation cannot be inte- 
grated analytically. Instead, we employ a shooting method to 
simultaneously solve Equations ( [26| , p7| ), and ([34]) subject 
to three boundary conditions: T = Tq and p = po at z = 0, and 

As in our non-radiative HBI simulations, we choose Q = 
2.0, Tz/To = 2.5, A) = 10^ and Kno' = 1500. Radiative cooling 
introduces another dimensionless free parameter. 



C 



~ Vth/H 



0.1 



0.01 cm- 



(35) 



We adopt a value of Co = 0.18, for which Wcooi — -0.1 Iwdyn at 
t = z = Q- Note that a simultaneous choice of Q, Kno, /3o, and 
Co implies specific dimensional values for our model cluster 
core: nj^o ~ 0.028 cm-^ kBTo ~ 3.0 keV, keTz ~ 7.5 keV, 
g ~ 1.5 X lO"** cm s-2, Bo ^ 0.26 pG, and Z ~ 248 kpc. ^ 

"* While this "central" density is a factor of ~2-3 smaller than those found 
in actual cluster cores, one cannot construct a thermodynamic equilibrium be- 
tween conduction and cooling with greater central densities (for our choices 
of f , Q, and Co). This is an indication that even unbridled conduction cannot 
offset radiative losses in all clusters (see, e.g., Zakamska & Narayan 2003 1. 




t= 50 




Figure 7. Magnetic-field lines (color-coded according to the local field 
strength) fro m r un H3dBrag at times / = 30 and 50 (in units of Hq/vh, q; 
see Equation |33| . The computational dom ain has dimensions L^xLy X L- = 
Ix 1x2 (in units of Hq; see Equation [32j. The magnetic-field strength is in 
units of (Avrpo)'''^. 



The corresponding units of length and time are [£] ~ 124 kpc 
and [f ] ~ 160 Myr, respectively; the former implies a grid size 
~0.24 kpc. The resulting thermodynamic equilibrium repre- 
senting our initial conditions is shown in Figure [8] in dimen- 
sional units. The temperature (solid line), ion density (dashed 
line), and inverse Knudsen number (dotted line) are very sim- 
ilar to those observed in the cool-core cluster A85, especially 
for r > 50 kpc (Cavagnolo et al.|2009| l. 
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Figure 8. Equilibrium atmosphere used as the initial condition in runs 
H2dBRad and H2dIRad. Profiles of the temperature (solid line), the ion den- 
sity (das hed line), and the inverse of the Knudsen number (dotted line; see 
EquationfTT) are given in the dimensional units on the accompanying scales. 
These prohles are very si milar to those observed i n the cool-core cluster A85, 
especially for r > 50 kpc |Cavagnolo et al.|2009) . 




10 20 30 40 
time 



10 20 30 40 
time 




20 30 
time 



Figure 9. Temporal evolution of the box-averaged horizontal (x) and vertical 
(z) kinetic and magnetic energy densities in runs H2dBRad (red lines) and 
H2dIRad (blue lines). The units of energy density and time ai'e, respectively, 
pov? „ ~ 2.7 X 10"'" ergs cm"^^ and ffo/vth.o - 160 Myr. 



Second, we must modify the temperature boundary condi- 
tion at z = in order to prevent the development of sharp tem- 
perature gradients between the ghost zones and the few first 
active zones where the cooling rate is greatest. We choose re- 
flective boundary conditions, which enforces a zero-gradient 
condition on the temperature at the bottom of the computa- 
tional domain. While this no longer implies a fixed heat flux 
through the computational domain (as in our non-radiative 
HBI runs), it is more compatible with the physical conditions 
in actual cluster cores. 

In Figure |9] we show the temporal evolution of the box- 



averaged kinetic and magnetic energy densities in runs 
H2dBRad (red lines) and H2dIRad (blue lines). In both runs, 
the instability growth rates are reduced from those in their 
respective non-radiative runs (see Figure [T]). The change 
in growth rates is due to an interplay between two effects. 
First, the equilibrium atmosphere has a shallower tempera- 
ture profile when cooling is taken into account. Because the 
HBI growth rate scales with (dlnT/dz)'/^, this naturally de- 
creases growth rates. On the other hand, as shown in Balbus] 
|& Reyn olds' (2008) by way of a local linear analysis of the ra- 
diative HBI, cooling acts to further destabilize the atmosphere 
(especially at longer wavelengths). Physically, this is because 
cooling weakens the ability of conduction to wipe away tem- 
perature fluctuations of a given wavelength along magnetic- 
field lines; such temperature fluctuations, a result of fluid el- 
ements' heat exchanges with the background conductive flux, 
are necessary for the HBI to function. The periodic oscilla- 
tions seen in the early evolution of the vertical kinetic energy 
are due to small departures from and oscillations about ther- 
modynamic equilibrium. As in our non-radiative HBI runs, 
Braginskii viscosity retards the development and growth of 
the HBI. The energy densities in runs H2dBRad and H2dIRad 
do not become comparable until t >4-Q (~6.4 Gyr). By the 
end of the simulations (f ~ 8 Gyr), the total magnetic energy 
in both runs has increased by a factor of ~'10. 

Figure [TO] exhibits the temperature (color) and magnetic- 
field linesXolack lines) in runs H2dBRad and H2dIRad at 
times f = 20, 25, 30, 40, and 50 (in units of //o/vth.o ^ 
160 Myr). As in the runs without cooling, the HBI devel- 
ops first at low altitudes and subsequently spreads to higher 
altitudes, with Braginskii viscosity significantly affecting the 
structure of the magnetic field for z > 0.2H(). This behavior is 
highlighted quantitatively in Figure 1 1 which shows the hor- 
izontally averaged magnetic -field angle {b^)x as a function of 
height at the same times as in Figure [Toj A comparison with 
Figure [3l reveals that radiative cooling further promotes hori- 
zontal augnment of the field lines. This is achieved not only 
by locally increasing the temperature gradient at the bottom 
of the box, which increases the local growth rate of the HBI, 
but also because the ensuing cooling flow helps to squeeze 
any horizontal field lines. Radiative cooling also introduces 
two other notable differences. 

First, thin (<1 kpc) filaments of cool (<1 keV) gas tran- 
siently appear throughout the course of run H2dBRad, typi- 
cally lasting anywhere between ^100 Myr and ^1 Gyr and 
extending over distances ~^30-80 kpc. Sometimes these fil- 
aments interact and merge with others nearby, forming rela- 
tively long chains of cool gas that wind throughout the cluster 
core. Sometimes filaments appear in pairs that run alongside 
one another over distances of a few tens of kpc. 

we focus in on three representative cool 
Figure 12 1 presents a w53 kpc x 74 kpc region 



In Figure 
filaments 

surrounding two neighboring cool filaments at time t 
7.6 Gyr and location z ~ 74—148 kpc. Figure \T2^ presents 
a w53 kpc X 53 kpc region surrounding a cool filament at time 
t « 8.7 Gyr and location z ~ 62-115 kpc. In all cases, these 
filaments follow the local magnetic-field lines, which serve to 
insulate them from the surrounding warm gas and which tend 
to become relatively isothermal over distances comparable to 
the Field length in the cool gas. This morphology is in agree- 
ment with dedicated numerical studies of local thermal insta- 
bility in globally stable, anisotropically conducting plasmas 
by |Sharma et al.| ( |20lO) l. The velocities in the filaments are 
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Figure 10. Spatial and temporal evolution of the radiative HBI with (top row) and without (bottom row) Braginskii viscosity. The temperature (color) and 
magnetic-field lines (black lines) are shown at times t = 20, 25, 30, 35, 40, and 50 (in units of //o/i'th.o — 160 Myr). The computational domain has dimensions 
Li-xLj = 1x2 (in units of Hg ~ 124 kpc). We have suppressed the imaging of temperatures beyond the fixed color-bar limits. 
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Figure 11. Spatial and temporal evolution of the horizontally averaged 
magnetic-field angle at each height in runs H2dBRad (red lines) and H2dIRad 
(blue lines). The units of length and time are, respectively, //o — 124 kpc and 
//o/i'th.o 160 Myr. 

also aligned with the local magnetic field and are well-ordered 
with speeds ^100-300 km s~', depending upon their orien- 
tation (vertically oriented filaments tend to have larger bulk 
velocities due to gravitational acceleration). The magnetic- 
field strength is enhanced during the formation of these fila- 
ments, reaching typical values of ^5-20 /iG. This enhance- 
ment extends over a region that is much longer than the ex- 
tent of the cold gas itself, since cool gas becomes compressed 
along field lines and evacuates regions of plasma. One con- 
sequence of this enhancement is the local production of hot 
(?s7-10 keV) gas due to parallel viscous heating, which often 
envelopes the cool filament and produces a sharper tempera- 



ture change across the filament. We discuss the astrophysical 
implications of these filaments in Section |7] as well as their 
agreement with current observational estimates. 

Second, while a cooling catastrophe inevitably occurs in 
both runs, the amount of time until the cooling catastrophe 
occurs, as well as the amount of cold mass at any given time, 
are different. In Figure [13] we show histograms of the frac- 
tions of the total mass in each thermal phase ('cool' refers 
to r < To ~ 3 keV and 'warm' refers to T > Tq) at times 
f = 20, 30, 40, and 50 (in units of //o/vth,o ^ 160 Myr) for 
runs M2dBRad (red) and M2dlRad (blue). Mass is binned by 
temperature in intervals of width A(T /To) = 0.1. The mass 
fraction in each phase is shown in each plot as a percentage. 
Braginskii viscosity delays a cooling catastrophe by reduc- 
ing the efficacy of the HBI. As a result, the amount of mass 
in the cool phase is always less when Braginskii viscosity is 
included. While the temperature distributions become com- 
parable by f = 50 (~8 Gyr), we note that in run M2dBrag (/) 
a significant portion of the cool gas is in the form of filaments 
located away from z = and (//) there is hot gas with temper- 
atures >7 keV due to parallel viscous heating. 

Why are there no cool filaments in run H2dlRad? The dif- 
ference has to do with the fact that local thermal instability 
does not gr ow exponentially in dy namically evolving atmo- 
spheres (e.g Balbus & Soker|1989 i. Without Braginskii vis- 
cosity suppressing the small-scale evolution of the HBI in the 
majority of the cluster core and thereby allowing conductive 
heating to offset a significant portion of the global radiative 
losses, a cooling flow readily develops and advects poten- 
tially unstable cold gas along with the bulk flow. Since a 
global equilibrium state is approximately preserved by potent 
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Figure 12. Representative cool filaments from run H2dBRad. (a) A 
Si53 kpcx74 kpc region surrounding two neighboring cool filaments at time 
t ^7.6 Gyr and location z ~ 74—148 kpc. (b) A ^53 kpc X 53 kpc region sur- 
rounding a cool filament at time t ^ 8.7 Gyr and location z ~ 62-1 15 kpc. In 
each of the frames, the temperature (color), magnetic-field lines (solid lines), 
and velocity vectors (black arrows) are shown; the temperature normalized 
to k^Tg ~ 3.0 keV. The magnetic field runs locally parallel to the filaments, 
with a strength ~5-20 fiG. The maximum magnetic-field strength and speed 
in each frame are (a) ~19.6 /iG and ~296 km s"', and (b) ~13.8 fiG and 
~250 km s"' , respectively. 

conduction outside of the innermost ^20% of the core, local 
thermal instability can proceed there. 

Finally, in Figure [14] we present the advective (red lines), 
convective (purple Imes), conductive (blue lines), and total 
(black solid lines) energy fluxes through z = 0.4 (~50 kpc) in 



runs H2dBR ad and H2dIRad. (See equations 20-22 of jBog- 
[danovic eF al. 2009 for definitions of these fluxes.) The black 
dashed lines denote the radiative cooling rate integrated over 
the region < z < 0.4, which must be balanced by the en- 
ergy fluxes for the core to be in equilibrium. Negative fluxes 
correspond to downward energy transport and vice versa for 
positive fluxes. Figure 14 shows that, in both cases, the atmo- 
sphere evolves from the mitial equilibrium at the time when 
the amount of inward conductive flux transported begins to 
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Figure 13. Fraction of the total mass in each thermal phase ('cool' refers 
to T < To ~ 3 keV; 'warm' refers to 7 > Tq) at times t = 20, 30, 40, and 
50 (in units of //o/v,h.o ^ 160 Myr) for runs H2dBRad (red) and H2dIRad 
(blue). The vertical dashed line denotes the temperature boundary separating 
the 'cool' and 'warm' phases. The total mass fraction in each of the phases is 
given as a percentage for each run. 
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Figure 14. Temporal evolution of the advective (red lines), convective (pur- 
ple lines), conductive (blue lines), and total (black solid lines) fluxes through 
z = 0.4 (~50 kpc) in runs H2dBRad and H2dIRad. The black dashed lines 
denote the radiative cooling rate integrated over the region < z < 0.4, which 
must be balanced by the energy fluxes for the core to be in equilibrium. Neg- 
ative (positive) fluxes correspond to downward (upward) energy transport. 
The units of time and energy flux are, respectively, ffo/^'th.o — 160 Myr and 



: 0.02 ergs s cm 



dwindle. The most striking difference between runs H2dIRad 
and H2dBRad is that the suppression of the conductive flux 
occurs much later in the run with Braginskii viscosity (t sa 20 
vs. f w 10). As a result, the conductive and total energy fluxes 
at the end of run H2dBRad are double those in run H2dIRad. 

After the onset of the HBI, run H2dBRad exhibits brief 
episodes of substantial convective flux, at times reaching 
^20% of the total flux. This reverse convective flux acts as 
a coolant in the energy equation. The convective flux max- 
ima correspond to the strong episodes of heat conduction to- 
wards the cool core, which can be traced to times when fila- 
ments cross the referent surface (z = 0.4) where the fluxes are 
evaluated. The presence of Braginskii viscosity promotes the 
formation of such filaments and, consequently, this mode of 
heat conduction and convective feedback. A possibility that a 
cool core mitigates abrupt changes in its thermal state via such 
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a feed back loop has been proposed by Balbus & Reynolds 
( 2008 1 as a mechanism to regulate thermal conductivity in hot 
(i.e. massive) galaxy clusters. 

6. MTI 

6.1. Background Equilibrium and Initial Perturbations 

We consider a non-radiative, plane-parallel plasma strat- 
ified in both density and temperature in the presence of a 
uniform gravitational acceleration in the vertical direction, 
g = -gz. We thread the plasma with a uniform background 
magnetic field oriented along x, so that there is no heat flux in 
the background state (i.e. q = 0). While a subcritical transition 
to turbulence exists for an initi ally vertica l magnetic field sub- 
ject to modest stirring (McCourt et al.|201 1 ), we focus only on 
the simplest background from which the linear MTI grows the 
fastest. Since q = Q, thermal equiUbrium is trivially satisfied 
and we are free to choose a linearly decreasing te mperature 
profile ( [Parrish & Stone|2005[ [McCourt et al.|20TT) i: 



T(z) = To 1 



Force balance then implies 



P(Z) = P()[ 1 



p(z) = Po 1 



3//n 



3Hn 



3Hq 



(36) 

(37) 
(38) 



Because the MTI induces large vertical displacements in the 
plasma, we atte mpt to minimize the effects of our bound- 
ary conditions (5 6.2 1 by sandwiching the unstable volume 
of plasma between two buoyantly neutral layers, following 
[Parrish & Stone] ( [2007[ ). We divide the vertical box size 
into three regions: region I (0 < z < L^/A-), region II 
(L,/4 <z< 3L,/4), and region III (3LJ4- <z< LJ. Re- 
gion II is described by Equations (36 1-(38 i with z — ?> z-L,/4, 
so that T = To, p = p^, and p = po at tEe base of the MTI- 
unstable region. Regions I and III are isothermal atmospheres 
with Ti = Tq and Tni = To (1 -L./6//0), respectively. Requiring 
force balance and continuity yields the following density and 
pressure distributions: 



Pi(z) 
Po 



Pi(z) 
po 



Pm(z) Pmiz) 



Po 



PO 



= Po 



= exp 
6H0 



Z-LJ4 
Ho 



exp 



Z-3LJ4 
Ho-Lje 



(39) 



(40) 



We further prescribe isotropic conduction in regions I and III 
in order to stabilize the MTI there. 

We apply Gaussian-random velocity perturbations to our 
background equilibrium, having a flat spatial power spectrum 
and a standard deviation of lO'^^v'th.o. While we have assumed 
equal ion and electron temperatures throughout this paper, we 
caution here that, due to the very low collisionality in cluster 
outskirts, the assumption of equal ion and electron tempera- 
tures may not hold. 

6.2. Numerical Setup and Boundary Conditions 

We non-dimensionalize our equations using the same units 
Note, however, that typical tem- 



chosen in Section 4.2 



peratures at the base of the MTI-unstable portion of the 



ICM, where the temperature begins to decrease outwards, are 
k^To « 5-8 keV, so that [i] - 400 kpc and [f] - 350 Myr are 
more representative numbers for the units of length and time, 
respectively. We choose the free parameters jio = 10^ and 
Kno' = 200, so that Bo ~ 0.016, v\\ ~ 2.4 x lO'^ T^l^p-\ and 

K|| ~ 0.12 T^l^p~^ in dimensionless units. While Kno under 
actual cluster-outskirt conditions is ^5 larger than our cho- 
sen value, we have found that the implied K|| (z = Lj, f = 0) ~ 1 
makes the computations unnecessaril y expensive. 

Although linear analyses without ( Balbus[2001[ l and with 
(Kll) Braginskii viscosity have shown that the fastest- 
growing MTI modes satisfy 1 ^ k\^H ^ /^'/^ and are there- 
fore l ocal, recent work on the nonlinear development of the 
MTI ( jMcCourt et al.[20rT| l has indicated that simulations with 
L/H ~ 0.1 significantly underestimate the magnitude of the 
turbulent velocities in the saturated state. We therefore choose 
box sizes HqxIHo (in 2d) and HoxHqxIHo (in 3d); the MTI- 
unstable layer then has a vertical size equal to Hq. 

Our boundary conditions are the same as in [McCourt et al.| 
( [201 1[ ). The temperatures at the upper and lower boundaries 
of our computational domain are fixed for all times, while 
the pressure is extrapolated into the upper and lower ghost 
zones in such a way as to ensure hydrostatic equilibrium at 
those boundaries. Periodic boundary conditions are imposed 
in the horizontal direction(s). We further maintain horizon- 
tal magnetic fields at the upper and lower boundaries on the 
computational domain. Note that this does not insulate the 
ghost zones from the computational domain, as we have pre- 
scribed isotropic conduction in regions I and III. Thus a fixed 
temperature difference is imposed across the vertical length 
of our model atmosphere. In this situation the MTI cannot 
exhaust the source of free energy, which is being constantly 
replenished by the boundary conditions and thus continuously 
drives MTI turbulence. If the temperature difference across 
the computation domain were allowed to relax (i.e. Neumann 
boundary conditions), the atmosphere woul d become isother- 
mal before the MTI could fully develop ([Parrish & Stone| 
[2005] [20071 [Parrish et al.|2()08i [McCourt et al.|201 1|l. 

Here we present results from five MTI simulations. 
M2dIsoP is a 2d simulation with isotropic pressure and serves 
as a reference run, enabling us to draw conclusions about the 
effects of Braginskii viscosity. M2dBrag is a 2d simulation 
with Braginskii viscosity and M2dBLim is a 2d simulation in 
which the pressure anisotropy is artificially limited using the 
Sharma et al. closure described in Section [T2] M3dBrag is 
a 3d simulation with Braginskii viscosity and run M3dIsoP 
is a 3d simulation with isotropic pressure. The parameters in 
these simulations are summarized in Table [T] 

6.3. 2d Simulations 

Figure [15] presents the evolution of the kinetic and mag- 
netic energies averaged over the MTI-unstable region (0.5 < 
z < 1.5). Runs in which pressure anisotropics are allowed 
to develop (red and purple lines) initially exhibit a growth 
rate equal to that of the run without Braginskii viscosity (blue 
line), in agreement with predictions from linear theory (Kl 1). 
This was expected on the grounds that the fastest-growing lin- 
ear MTI modes are Alfvenically polarized (i.e. = 0) and 
therefore escape viscous damping by not producing a linear 
pressure anisotropy. However, these modes do produce a non- 
linear pressure anisotropy, which begins to significantly de- 
crease the growth rate once (SBj_/B f > 0.1. The growth rate 
in run M2dBLim (purple fine) does not show this behavior. 
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Figure 15. Temporal evolution of the horizontal (x) and vertical (z) kinetic 
and magnetic energies, averaged over the iVITI-unstable portion of the box 
(0.5 <z< 1.5), in runs IVIldBrag (red lines), M2dIsoP (blue lines), and 
M2dBLim (purple lines). The units of e nerg y density and time are, respec- 
tively, Po^i ^o/vth.o (see EquationF 
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Figure 16. Temporal evolution of the magnetic-field angle i)?, averaged over 
the MTI-unstable portion of the box (0.5 < z < 1.5), in runs M2dBrag (red 
line), M2dIsoP (blue l ine) , and M2dBLim (purple line). The unit of time is 
^^0 / I'th.o (see EquationL 

instead agreeing with that of run M2dBrag (red line). This is 
because the limiters do not allow the local nonlinear pressure 
anisotropy to become greater than the local magnetic pres- 
sure. All three runs reach an approximately saturated kinetic 
energy corresponding to a Mach number of several percent. 
The total magnetic energy increases by a factor of ^100 over 
the course of the runs, with a final /? ~ 10^ and most of the 
energy in the vertical component. The kinetic and magnetic 
energies saturate in approximate equipartition. 
Curiously, runs in which the pressure may become 



anisotropic have a slightly prolonged phase of exponential 
growth and, consequently, greater maximum energy densi- 
ties than those found in run M2dIsoP. We believe there are 
two principal reasons for this difference. First, Braginskii vis- 
cosity suppresses the formation of small perpendicular scales, 
which results in conditions less favorable for grid-scale mag- 
netic reconnection. Second, the suppression of perpendicu- 
lar fluctuations due to Braginskii viscosity constrains the per- 
turbed magnetic field to grow predominantly vertically in run 
M2dBrag. As a result, there is less interference amongst the 
buoyant plumes than in run M2dIsoP, where the perturbed 
field lines acquire a more tangled topology. Indeed, runs 
M2dBrag and M2dBLim exhibit greater saturated values of 
/^^^^ than those in run M2dIsoP (see Figure 
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These differences are dramatically highlignted in Figure 17 
which shows the overall spatial and temporal evolution of 
the atmosphere during each of our 2d runs. The temperature 
(color) and magnetic-field lines (solid lines) are displayed in 
each of the six frames, which show the atmosphere at the dif- 
ferent times f = 17.5, 20.0, 22.5, 25.0, 27.5, and 30.0 (in units 
of //()/vth,o). Clear differences exist in the topology and evo- 
lution of the magnetic field between each of the runs. 

The magnetic -field fluctuations in runs M2dBrag and 
M2dBLim emerge on larger scales than in run M2dIsoP. This 
is because the initial growth of modes with (5B|| ^ induces a 
positive pressure anisotropy that shifts all unstable modes to 
larger (parallel) wavelengths. (Recall from Equation 25 that, 
for our chosen plasma parameters, fluctuations in magnetic- 
field strength as small as ^0.1% can produce a pressure 
anisotropy comparable to the magnetic tension.) In addi- 
tion, the growing modes in runs M2dBrag and M2dBLim ex- 
hibit more of a sawtooth-like structure than do those in run 
M2dIsoP. This is because Braginskii viscosity targets only 
those motions that change the field strength. As a result, 
magnetic-field lines tend to remain locally parallel to one an- 
other in order to minimize field-line compressions and rar- 
efactions. This also results in a more laminar development 
of the instability than in run M2dIsoP, as small-scale features 
along field lines are viscously damped. 

Whether the pressure anisotropy is limited (M2dBLim) or 
not (M2dBrag), the turbulence tends to arrange the magnetic 
fields in long, thin flux sheets with ^ the field reverses 
its direction at the smallest scale available to it (the numer- 
ical resistive scale) but field lines curve at the scale of the 
flow (the viscous scale) except in the sharp bends o f the folds 
(for an analytical theory of folded structure, see Schekoch"P] 
|hin et al. 20 02 ). This folded structure, a general property or 
random forcing in plasmas with large magnetic Prandtl num- 
bers, allows the small-scale direction-reversing magnetic field 
to back react on the flow in a spatially coherent way: the ve- 
locity gradients become locally anisotropic with respect to 
the direction of the folds (with bb.'Vv partially suppressed) 
so that the dynamo saturates at marginally stable balance be- 
tween reduced "parallel" stretching and "perpendicular" mix- 
ing. Thus, the fie ld strength and the field-l ine curvature are 
anticorrelated (see'Schekochihin et al.|2004| . 

One consequence of reduced parallel stretching is the no- 
table paucity of firehose instabilities in run M2dBrag. In 
run H2dBrag (HBI with Braginskii viscosity), such a reduc- 
tion in parallel stretching and concomitant production of mi- 
croscale instabilities was not possible, since correlations be- 
tween field strength and field-line curvature are necessary for 
the HBI to function in the first place. Rather than reduce the 
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Figure 17. Spatial and temporal evolution of the MTI with Braginskii viscosity (top row), without Braginskii viscosity (middle row), and with limited Braginskii 
viscosity (bottom row). The temperature (color) and magnetic-field lines (black lines) are sho wn at times t = 17.5, 20.0, 22.5, 25.0, 27.5, and 30.0; see Equation 
^3) . The computational domain has dimensions L^xL^ = 1x2 (in units oi Hq; see Equation[32j. We have suppressed the imaging of temperatures beyond the 
nxed color-bar limits. 



pressure anisotropy to marginally stable values (Equation 24 1 
via the secular growth of microscale fluctuations, MTI-driven 
turbulence appears to avoid the production of large pressure 
anisotropies altogether by orienting magnetic fields primar- 
ily across local velocity gradients. Consequently, the box- 
averaged pressure anisotropy (not shown) is sub-marginal 
during the nonlinear phase of evolution. 

6.4. 3d Simulations 

Figure [18] shows the evolution of the kinetic and magnetic 
energies averaged over the MTl-unstable region (0.5 < z < 
1.5) for runs M3dBrag (red lines) and M3dlsoP (blue lines). 
Many features are in agreement with our 2d simulations. First, 
the growth rates from both runs are equal until the nonlin- 
ear pressure anisotropy begins to significantly decrease the 
M3dBrag growth rate. Second, run MBdBrag has a slightly 



prolonged phase of exponential growth and greater maximum 
energy densities than those found in run M3dIsoP. As ex- 
plained in Section |6.3[ we believe this is due to Braginskii 
viscosity suppressing the formation of small perpendicular 
scales. While a certain amount of the decrease in saturated en- 
ergies from 2d to 3d is likely due to the decreased resolution, 
it is also due to the fact that, in 3d, magnetic field lines are al- 
lowed to penetrate into the additional dimension and alleviate 
regions of strong magnetic pressure. Overall, the evolution of 
the average energies are qualitatively similar in 2d and 3d. 

The average magnetic-field angle, on the other hand, shows 
qualitative differences between runs M3dBrag and M3dlsoP 
(Figure 19 1. In run M3dBrag the magnetic field at the end 
of the Imear phase is oriented slightly more horizontal than 
in run M3dlsoP. However, subsequently grows faster and 
ultimately becomes larger than in run M3dIsoP. We attribute 
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Figure 18. Temporal evolution of the horizontal ("h") and vertical ("z") ki- 
netic and magnetic energy densities, averaged over the MTI-unstable portion 
of the box (0.5 < z < 1.5), in runs M3dBrag (red lines) and M3dIsoP (blue 
lines). The units of energy density and time are, respectively, poifdo 
^^o/i'th.o (see Equation|33|. 
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Figure 19. Temporal evolution of the magnetic-field angle fo?, averaged over 
the MTI-unstable portion of the box (0.5 < z < 1.5), in runs M3dBrag (red 
line) and M3dIsoP (blue line). The unit of time is Hq / Vth,o (see Equation|33|. 



this difference to the presence of the Alfvenic MTI: as the 
mean field becomes more vertical, Braginskii viscosity cou- 
ples ky y Alfven modes to damped slow modes and drives 
them buoyantly unstable. These modes, absent in the case of 
isotropic pressure, likely play a role in reorienting the mag- 
netic field vertically at late times. 

Figure 20 exhibits the magnetic-field strength (color) at 
time t = 37 in runs M3dBrag and M3dIsoP. As in runs 
M2dBrag and M2dIsoP, the magnetic-field fluctuations in run 
M3dBrag emerge on larger scales than in run M3dIsoP. In ad- 




MSdlsoP 




Figure 20. Pseudocolor plot of magnetic-field strength (color) in runs 
M3dBrag (ton) and M3dIsoP (bottom) at time t = 37 (in units of //o/»'th.o; 
see Equation |33| . The computational dom ain has dimensions L^xLyXL- = 
Ix 1x2 (in units of Hq; see Equation [32j. The magnetic-field strength is in 
units of (4npo)^/^. 



dition, because small-scale features along field lines are vis- 
cously damped, there is less reconnection and the magnetic 
field remains coherent over longer distances. As a result, the 
folded structure of the field is more apparent in run M3dBrag. 
It is also clear from this figure that the field strength and the 
field-line curvature are anticorrelated. 

7. DISCUSSION 
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7.1. HBI 

Our simulations of the HBI that self-consistently allow for 
anisotropic pressure have a number of implications for the 
structure and evolution of cool-core clusters, some of which 
may be potentially observable. Perhaps the most important 
is our finding that field-line insulation of the entire cool core, 
previ ously found to be a nonlinear consequenc e of the HBI 
(e.g. Parrish et al. 2009 ; Bogdanovic et al. 2009 ), does not oc- 
cur on astrophysically relevant timescales. The smaller degree 
of collisionality outside of the innermost regions of cool-core 
clusters ensures that the magnetic-field lines retain a strong 
vertical component outwards of ^50 kpc from the cluster cen- 
ter There, the comparatively large pressure anisotropy self- 
generated by the HBI forces the fastest-growing modes to 
have relatively long wavelengths, at w hich non-local effect s 
start to play a role and curb growth (see Latter & Kunz 2012| l. 
However, the relatively low temperatures and high densities 
in the innermost regions result in very little difference there 
between the inviscid and viscous cases. In the presence of 
radiative cooling, a cooling catastrophe inevitably occurs in 
these regions. This highlights the need for radio-mode feed- 
back at these scales from a powerful central dominant galaxy. 

Another important result from our simulations is the for- 
mation of cool filaments when anisotropic conduction, Bra- 
ginskii viscosity, and radiative cooling are all taken into ac- 
count. The pressure anisotropics generated by the HBI sup- 
press its ability to impede the conductive flux to small radii. 
As a result, much of the atmosphere evolves slow enough to 
allow local thermal instability to set in. The physical char- 
acteristics of the magnetically aligned, cool filaments that 
emerge in our simulation (temperatures <1 keV; magnetic- 
field strengths ^5-20 /iG; magnetically aligned velocities 
--100-300 kms"'; hfetimes ~10^-10^ Myr; aspect ratios 
^30-80) ar e similar to those observed or observationa lly in- 
ferred (e.g. [Hatch et al.|[2006} [McDonald eTd?|[20T^ . The 
magnetic field is responsible for insulating the filaments from 
the surrounding hot gas, with strengths capable of stabilizing 
the filaments and possibly delaying star formation (Fabian et^ 
[^[2008| l. In other words, our simulations are able to repro- 
duce many of the observed properties of the cool filaments 
without a need for ad hoc heating presc riptions to offset ra- 
diative cooling (such as those used by McCourt et al.[[20T? 



and Sharma et al. 2012). However, we caution that our simu- 
lations may be overestimating the thicknesses of the filaments 
since (/) we have neglected line cooling below T ^ 4 x 10^ K 
and (//) our spatial resolution is only ~0.24 kpc. As a result, 
we are not able to form the very thin (<100 pc) filaments that 
are sometimes observed in, for example, NGC 1275 (Fabian 
[et al.|20 08 ) and Abell 1795 ( McDonald et al.|20lO| ). 

We have also found that the cool filaments formed in our 
simulations are often surrounded by gas that is hotter than av- 
erage. This is because the compression of field lines during 
the formation of the filament induces a pressure anisotropy 
that leads to parallel viscous heating as it is collisionally re- 
laxed. Temperatures in these hot filamentary "envelopes" are 
in the range «7-10 keV. It would be interesting to examine 
deep Chandra observations of nearby cool-core clusters for 
evidence of hot filaments. Since these hot filaments would be 
relatively tenuous and hence have low emissivity, they could 
easily have been missed in existing analyses. However, the 
hot filaments may be revealed by constructing maps of the 
Fe XXVI/Fe XXV K-shell recombination line ratio. 

Finally, the structure of the magnetic field produced by the 



HBI is significantly affected by Braginskii viscosity. Bragin- 
skii viscosity causes the HBI to grow efficiently only for a 
thin band of (parallel) wavelengths, which correlate with the 
local collisionality of the plasma (recall that modes with non- 
negligible growth rates satisfy k\\H < \/Kn < 3k\\H\ Kll). 
The innermost (relatively collisional) regions of cool-core 
clusters are thus likely to harbor preferentially azimuthal mag- 
netic fields due to efficient field-line reorientation by the HBI. 
However, AGN- and/or merger-driven turbu lence may be able 
to randomiz e this field (Parrish et al. '20l0; Ruszkowski & Ohj 
,2010, McCourt et al.|2011| ). At distances -50-100 kpc from 
the cluster center, we have found that the flow of heat oc- 
curs primarily along magnetic sheaths (or filaments). Beyond 
— 100 kpc (where the collisionality is low) our results suggest 
that the HBI exerts relatively little influence on the direction 
of the magnetic field over astrophysically relevant timescales. 
These findings can be further tested by well-resolved, global 
simulations that include anisotropic conduction, Braginskii 
viscosity, and turbulent stirring, as well as by cluster obser- 
vations with the Expanded Very Large Array and (eventually) 
the Square Kilometer Array using rotation measures of back- 
ground polarized radio sources (e.g. jBogdanovic et al.|201T| ). 

7.2. MTI 

The consequences of including Braginskii viscosity in a 
treatment of the MTI are not as great as for the HBI. As 
a result, many of the observationally important results pre- 
viously found in MTI simulations that neglected Braginskii 
viscosity (e.g. radially biased magnetic fields, flattened tem- 
perature profiles, vigorous subsonic turbule nce that may lend 
an appreciable a mount of pressure support; 'Parri sh & Stone 
2005 2007 ; Parr ish et al.( 2008 ; McCourt et al. 201 Ij [Parrish 
|et al.||201 1| ) are for the most part unchanged. The reason is 
simple: the fastest-growing MTI modes are polarized such 
that 5B\\ =0 and therefore escape parallel viscous damping. 
There are some differences, however, that concern the emerg- 
ing structure of the magnetic field. 

For example, the magnetic-field fluctuations produced by 
the MTI with Braginskii viscosity emerge on larger scales 
than do those without Braginskii viscosity. This is because the 
generation of a positive pressure anisotropy shifts all unstable 
modes to larger (parallel) wavelengths. In addition, because 
field-line compressions and rarefactions are damped by the 
linear pressure anisotropy they generate, magnetic-field lines 
tend to remain locally parallel to one another. As a result, the 
magnetic field remains coherent over larger distances. Since 
the viscosity of the ICM is much larger than the magnetic re- 
sistivity, the turbulence tends to arrange the magnetic fields 
in long, thin flux sheets with <C i\\ — ^visc, for which the 
field strength and the field-line curvature are anticorrelated 
(see Schekochihin et al. 2004 ). 

It follows that Braginskii viscosity implies a minimum 
physical scale above which the magnetic field can fluctuate 
freely. This is interesting in light of observations of the Fara- 
day rotation measure, which indicate magnetic-field correla- 
tion lengths in the range —0. 1 to a few tens of kpc (e.g. Vogt & 
'EnBHn 2005 ; Guidetti et al. 2008 ; Bonaf ede et al.|20101[Gov- 
,oni et al. 2010, Kuchar & EnBlin 2011:> . While Braginskii 
viscosity may not be unique in its effect on the geometry of 
the intracluster magnetic field, it may provide important clues 
about the physical processes controlling its correlation length. 
On the other hand, if Braginskii viscosity can be shown to be 
the dominant mechanism, it would provide a direct link be- 
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tween the macroscopic observables such as the field correla- 
tion length and the elusive, microscopic plasma processes. 

7.3. Comparison with Related Work 

Shear viscosity, whether isotropic or anisotropic, has long 
been suspected to play an important role in the ICM. Using 
Chandra X-ray observations of Ha filaments in the Perseus 
cluster, Fabian et al. (2003b) argued that the effective viscos- 
ity of the ICM must be large enough to explain the appar- 
ently laminar flow as rising bubbles drag up colder, inner gas. 
[Fabian et al. ( |2003aj l also showed that viscosity can have an 
important effect in dissipating the sound energy produced by 
the formation of bubbles and heating the surrounding ICM. 
Subsequent numerical work by Reynolds et al. (2005 ) demon- 
strated that viscosity may be necessary to mamtam the ob- 
served integrity of AGN -blown buoyant cavities by quenching 
Kelvin-Helmholtz and Rayleigh-Taylor instabilities. Elimi- 
nating the need for bubbles to inflate supersonically in order 
to evade these instabilities, viscosity also provides a natural 
explanation for the obse rved absence of strong shocks bound- 
ing the ghost cavities. Dong & Stone ( 2009 1 extended this 
work by taking into consideration Braginskii viscosity and 
studyin g the influence of d ifferent magnetic-field orientations. 
Finally, [Kunz e t ar](|2011) showed that Braginskii viscosity, 
when regulated by microscale instabilities, provides a local 
thermally stable heating source. Given a sufficient supply of 
turbulent power, this provides a physical mechanism for miti- 
gating cooling flows and preventing cluster core collapse. Our 
study compliments all these efforts and lends credence to the 
notion that understanding the viscosity of the ICM is vital to 
understanding its morphology, energetics, and stability. 



In parallel to the work presented here, Pa rrish et al.| ( 2012] 
hereafter P12) carried out an independent study of the HB 
and MTI subject to Braginskii viscosity. In the areas of over- 
lap between our work and that presented in P12, there is broad 
agreement. However, there are some differences worth not- 
ing, which we believe are mainly due to different choices of 
free parameters, initial conditions, and numerical approaches. 

As discussed in Section 3.3 our choice of free parameters 
is motivated by the physical conditions in galaxy clusters and 
the numerical constraints related to the impact of microscale 
instabilities on our results. Accordingly, we have chosen to 
initialize the magnetic-field strength in all of our simulations 
using /3 ^ lO'^-lO^. In all cases, the energy in the magnetic 
field saturates with /? ~ 10^-10^, in approximate equiparti- 
tion with the fluid motions. By comparison, in P12 the local 
HBI and MTI simulations have an initial (3 ^ 10'^, the fidu- 
cial global HBI simulations have an initial /3 ^ 10^, and the 
global MTI simulations have an initial /3 ~ 10^-10^. None 
of the runs presented in P12 appear to saturate with approx- 
imate equipartition between kinetic and magnetic energies. 
One consequence of the different choice of /? is that the satu- 
rated magnetic field in our MTI simulations remains strongly 
biased in the vertical direction, whereas in the P12 simula- 
tions the magne tic field be comes nearly isotropic (a result also 
found by McCo urt et al.| [201 1 ). In the saturated state of our 
MTI runs, the horizontal kinetic energy is less than the verti- 
cal magnetic energy and so there is insufficient energy in the 
horizontal motions to isotropize the magnetic field. 

Differences between our results and those in P12, especially 
concerning the evolution of the HBI, may also be due to dif- 
ferent choices of Pr. We have used Pr ~ 0.02 (see Equation 
fTSll, while P12 chose Pr = 0.01. Our non-radiative and ra- 



diative HBI simulations also start with atmospheres (similar 
to A 1795 and A85) that are less coUisional than the fidu- 
cial cool-core model in P12. In the notation of P12's figure 
2, our cool-core atmospheres have Wcond/'^buoy ^ 10-30 for 
r > 20 kpc (i.e. Wvisc/^buoy ^ 2-5). This may be the reason 
why we find the HBI to be appreciably suppressed beyond 
^50 kpc, while those authors do not. However, this may also 
be due to the fact that all of our HBI simulations started with 
a magnetic field aligned with gravity, whereas P12 initialized 
their simulations with tangled magnetic fields on scales 30- 
50 kpc. Future work may resolve these differences. 

Finally, in our paper we have highlighted the influence of 
Braginskii viscosity on the morphology of the intracluster 
magnetic field and on the formation of cool filaments. We 
have also tried to address the impact of microscale instabilities 
on our results by taking two different approaches to capture 
their macroscale effects: (/) by working at very high spatial 
resolution so that the microscale instabilities grow fast enough 
to naturally regulate the pressure anisotropy, and (//) by em- 
ploying anisotropy limiters to restrict the pressure anisotropy 
to stable or marginally stable values. Using these approaches, 
we have shown that the manner in which microscale instabili- 
ties saturate affects the properties of the intracluster magnetic 
field. These approaches must be considered provisional, how- 
ever, as there is currently no complete microphysical theory 
concerning the saturation of these instabilities. 

7.4. Summary and Outlook 

In this paper, we have employed numerical simulations to 
investigate the linear and nonlinear dynamic and radiative 
stability of a weakly collisional, magnetized ICM. We have 
taken into consideration the effects of anisotropic heat and 
momentum transport, radiative cooling, magnetic tension, and 
microscale instabilities, and have ascertained a number of 
their implications for the structure of the intracluster magnetic 
field, the resolution of the cooling-flow problem, and the na- 
ture of convective turbulence in a dilute plasma. 

Despite such progress, there are still a number of unan- 
swered questions, some of which may be addressed by well- 
resolved global, three-dimensional numerical simulations of 
cluster cool cores and cluster outskirts. However, the efficacy 
of such simulations is likely to be contingent upon the im- 
plementation of a realistic sub-grid model for the microscale 
instabilities that captures their interplay with the computation- 
ally resolved meso- and macroscales. While formulating such 
a model is a rather formidable task, dedicated efforts to con- 
struct a more complete microphysical theory and to under- 
stand its bearing on heat and momentum transport, magneto- 
genesis, and thermodynamic stability in astrophysical systems 
are clearly needed. 
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